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ABSTRACT 


Daily  observations  of  sea-surface  temperature  (SST)  taken  just  north 
of  Pt.  Sur,  California  over  the  past  12  years  are  examined  for  statistical  and 
oceanographic  content.  Spectral  analysis  of  the  data  reveals  an  annual  cycle 
and  two  harmonics,  but  no  obvious  influence  from  the  local  winds  or  tides.  A 
statistical  model  has  been  constructed  based  on  a  decomposition  of  the  data 
into  its  various  deterministic  and  random  components.  As  expected,  the  model 
provides  forecasts  which  are  significantly  better  than  those  that  would  be 
obtained  from  a  simple  cl imatological  treatment  of  the  data  alone.  The  model 
is  also  used  to  generate  simulated  series,  which  in  turn  are  used  as  a  basis 
to  estimate  adequate  sampling  rates. 

In  6  years  of  the  12  year  record,  an  abrupt  major  decrease  in  SST  was 
observed  between  February  and  April  and  has  been  identified  as  the  spring 
transition  to  coastal  upwelling.  One  such  transition  was  found  to  occur 
almost  simultaneously  at  three  locations,  including  Pt.  Sur,  along  the 
California  coast  between  Pt.  Conception  and  San  Francisco.  The  data  also 
indicate  that  a  well-defined  annual  temperature  cycle  does  exist  at  Pt.  Sur 
although  its  range  is  reduced  and  its  maximum  occurs  later  in  the  year  than 
would  normally  be  expected  at  this  latitude.  Three  El  Nino  events  are  clearly 
evident  in  the  12-year  record  including  the  present  El  Nino;  SSTs  have 
increased  by  almost  3C  at  Pt.Sur  as  a  result  of  the  present  El  Nino  warming. 
Over  the  whole  12  year  period  an  average  increase  in  temperature  of  1.6C  was 
found. 
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ANALYSIS  OF  A  12-YEAR  RECORD  OF  SEA-SURFACE 
TEMPERATURES  OFF  PT.  SUR,  CALIFORNIA 


I.        Introduction 

Sea-surface  temperatures    (SSTs)  are  collected  on  a  regular  basis    (SIO 
Reference  81-30)  along  the  California  coast  at  approximately  2b  locations.      In 
view  of  the  increasing  interest   in  coastal   circulation  with   respect  to  the 
dispersal   of  pollutants  and  a  continuing  interest  in  fisheries   biology,    it   is 
surprising  that  few  systematic  studies  have  been  conducted  to  examine  these 
records,   particularly  within  the  past  15  years.     Perhaps  the  availability  of 
high-quality  temperature  data  along  the  continental   shelves  from  new  and 
sophisticated  instruments  has  caused  us  to  overlook  this   resource  in  recent 
years.     However,   the  coastal   observations,   unlike  those  acquired  from  moored 
arrays,   often  extend  over  many  years  and  thus  provide  a  unique  opportunity  to 
examine  coastal    variability  over  relatively   long  periods.     For  locations  where 
the  measuring  site  has  a  good  exposure  to  the  adjacent  continental   shelf  and 
slope,   measurements  of  SST  may  be  particularly  revealing  with   respect  to  some 
of  the  physical   processes  that  occur  regionally  as  well   as   locally.      In  this 
regard  a  12-year  record  of  SSTs  off  the  central   California  coast  at  Granite 
Canyon,    15km  south  of  Monterey,   has  been  chosen  for  the  initial   study.     These 
data  have  been  acquired  daily  by  the  Marine  Culture  Laboratory  of  the 
California  Fish  and  Game  Commission  since  March  of  1971.     The  Granite  Canyon 
site  where  the  observations  are  made   (see  Fig.    1)  has   an  excellent  exposure  to 
the  deep  ocean  with  the  continental   shelf  extending  less  than  10km  offshore. 

Previous   investigations  along  the  U.S.  West  Coast  indicate  that  the 
oceanographic  information  content  of  coastal   SST's   is  usually  high. 
SSTs  have  been  shown     (a)  to  be  useful   indicators  of  ocean  temperature 
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FIGURE      1 


Measurement  Locations.      The  principal   time  series   considered  here 
was   taken  at  Granite  Canyon.     Qualitative  comparisons   are  made  to 
fragmentary  series   of  subsurface  temperatures   close  by,    and  to  a 
short  time  series  at  Diablo  Canyon.     Comparisons   of  the  Granite 
Canyon  SSTs  will    be  made  to  extensive  time  series   of  SSTs   at   the 
Hopkins  Marine  Station   at   Pacific  Grove  and  to  SSTs   at   the 
Farallon   Islands. 
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variability,  (b)  to  be  often  representative  of  phenomena  occurriny  over  wide 
regions,  (c)  to  be  related  to  other  ocean  and  atmospheric  variables,  and  (d) 
to  have  consistent  internal  structures.  Based  on  surface  and  subsurface 
temperatures  along  the  California  coast,  Reid,  Roden  and  Wyllie  (1958) 
indicated  two  different  causes  of  seasonal  variability  in  SSTs.  The  first  is 
coastal  upwelling.  North  of  Pt.  Conception  the  upwelling  process  decreases 
the  seasonal  range  in  temperature  and  increases  the  duration  of  cooler 
temperatures.  The  second  cause  stems  from  the  importance  of  the  coastal 
countercurrent.  This  flow  to  the  northwest  transports  warm  water  from  the 
south  along  the  California  coast,  increasing  surface  temperatures  during  the 
fall.  Also,  Reid,  Roden  and  Wyllie  (1958)  found  that  variations  in 
temperature  along  the  U.S.  West  Coast  show  considerable  coherence  over  great 
distances  and  that  major  anomalies  in  temperature  persist  for  several  months. 

Using  coastal  SSTs  (and  sea  level)  Stewart*  Zetler,  and  Taylor  (1958) 
verified  the  existence  of  anomalous  conditions  off  the  U.S.  West  Coast  during 
1957  and  1958.  Robinson  (1957)  indicated  the  value  and  representativeness  of 
coastal  SSTs.  She  found  that  coastal  SSTs  and  air  temperatures  were  highly 
correlated  and  that  SST  anomalies  identified  at  coastal  locations  could  be 
correlated  with  SST  anomalies  up  to  250km  offshore.  Roden  (1963)  found  that 
non-annual  extreme  temperature  fluctuations  were  coherent  and  in  phase  over 
distances  of  200-30Ukm.  He  also  found  that  the  average  persistence  of 
non-annual  minimum  temperatures  was  3-4  months  while  the  average  persistence 
of  non-annual  maximum  temperatures  was  3-5  months  along  the  coast  of  northern 
California.  Spectral  analysis  indicated  the  presence  of  an  annual  cycle  only, 
however. 

More  recently  List  and  Koh  (1976)  analyzed  SSTs  along  the  California 
coast  over  several  years  using  spectral  analysis  and  found  that  many  of  the 


coastal   stations   exhibited  strong  spectral   peaks   at   14.7  days   and  that  these 
peaks  were  probably   the  result  of  tidal   aliasing  and/or  the  spring  and  neap 
tides.     They  also  found  that   a  well-defined  seasonal    variation   in  SST  was 
essentially  absent   between  Pt.   Conception  and  Monterey  Bay  and  that   certain 
fluctuations  occurred  simultaneously  at  most  coastal   stations. 

Surface  and  subsurface  temperatures   across  Monterey  Bay  were  collected  on 
a   regular  basis   by  the  Hopkins  Marine  Station  from  1929  to   1974.     Results   of 
these  measurements   have  been   reported  by  Skogsberg   (1936),   Skogsberg  and 
Phelps    (1946),    and  Bolin  and  Abbott    (1963).     Much  of  the  data  presented  in 
these  references  were  time-averaged  providing  monthly  mean  thermal   conditions 
over  the  Bay.     Skogsberg  used  these  data  to  define  three  separate  oceanic 
seasons  for  the  Monterey  Bay  area. 

Our  purpose  in  this  paper  is  threefold:      first,   to  model    (i.e.   describe 
mathematically  and  statistically)  the  behavior  of  the  12  year  Granite  Canyon 
data  set;   second,   to  use  the  model   and  other  statistical   techniques  to  project 
or  predict  the  data  and  estimate  time  scales  of  variability  in  SSTs,   and 
third,   to  provide  a  descriptive  interpretation  of  the  Granite  Canyon  data  from 
an  oceanographic  viewpoint. 

2.       The  Granite  Canyon  Data 

SSTs  at  Granite  Canyon  are  measured  daily  at  approximately  0800 
local   time.      The   raw  data  are  shown  in  Figure  2.   Because  the  data  are  taken 
only  once  every   24  hours,   the  possibility  of  aliasing  by   the  predominant 
semidiurnal   tide  exists.     This  problem  will   be  addressed  later  in  the 
analysis.     As  mentioned  before  the  location  has   an  excellent   exposure  to  the 
deep  ocean.     Temperatures  are  read  to  the  nearest  0.1C  using  a  calibrated  VWR 
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thermometer,  and  are  considered  accurate  to  about  ±  0.2C  (SIO  Reference  81-30). 
The  observations  span  12  years  starting  March  1,  1971  and  ending  on  February 
28,  1983.  SSTs  over  the  same  time  period  were  acquired  for  the  Hopkins  Marine 
Station  at  Pacific  Grove  and  the  Farallon  Islands  (see  Figure  1  for  location). 
The  accuracy  of  these  data  is  also  about  ±  0.2C. 

Since  the  data  from  Granite  Canyon  (GC)  have  only  been  available  since 
1971  and  since  they  have  not  been  compared  with  other  coastal  temperatures,  it 
was  of  interest  initially  to  consider  the  representativeness  of  these 
observations.  Consequently  we  have  compared  the  Granite  Canyon  data  with 
several  other  data  sets  for  periods  where  comparisons  could  be  made.  The  annual 
variations  in  SST  were  generally  similar  at  Granite  Canyon,  Diablo  Canyon  and 
the  Farallons,  but  temperatures  were  consistently  higher  at  Pacific  Grove 
during  the  spring  and  summer  months  than  they  were  at  the  other  three 
locations. 

Temperature  data  at  four  depths  between  113  and  510m  were  acquired  from  a 
moored  array  approximately  50km  south  of  Granite  Canyon  for  the  period  03/05/80 
to  04/07/80  and  compared  with  the  SST  data  at  Granite  Canyon,  see  Figure  3. 
The  period  was  chosen  because  it  includes  a  major  event  tentatively  identified 
as  the  spring  transition  to  upwelling,  a  phenomenon  which  has  been  recently 
identified  further  north,  along  the  coast  of  Oregon,  by  Huyer  et  al,  (1979). 
This  event  can  be  readily  identified  in  a  number  of  the  annual  cycles  in  the 
raw  data  (Figure  2,  and  Table  1).  The  drop  in  temperature  associated  with  this 
event  during  March  of  1980  can  also  be  seen  at  each  depth  in  Figure  3,  down  to 
510m. 

Figures  4,  5,  6  show  comparisons  of  SSTs  between  Granite  Canyon  and 
Diablo  Canyon,  between  Granite  Canyon  and  Pacific  Grove,  and  between  Granite 
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Canyon  and  the  Farallon   Islands   respectively   for  the  year   1980.     Diablo  Canyon 
is    located  about  200  Km  south  of  Granite  Canyon  and  the  data  from  that   location 
were   recorded  at   5m  below  the  surface  on  the  continental   shelf,   just  off   the 
PG&E  nuclear  power  plant   located  on  the  adjacent  coast.     The  temperatures  at 
Pacific  Grove  are  collected  near  the  shore  in  a  sheltered  area  inside  Monterey 
Bay.     The  Farallon   Islands   data  contained  numerous  multi-day   gaps  which  were 
filled  by   linear  interpolation.     The  same  spring  transition  shown  in  Figure   3 
is  seen  to  be  present   in  these  figures.     The  overall   patterns   in  the  annual 
cycle  at  these  four  locations  are  seen  to  be  generally  similar  at   least  during 
1980. 

3.       Gross   trends   in  the  Granite  Canyon  SST's 

We  discuss   here  gross   features   of  the  Granite  Canyon  SST   record.     These 
features  are  of  interest  per  se  and  are  also  a  prelude  to  the  more 
quantitative  decomposition  and  time  series  analysis  performed  in  the 
subsequent  sections. 

Figure   7  shows   a  scatter  plot  of  the  temperature  data  against  day   number. 
In  most  years  fewer  than  three  observations  were  missing  from  the  raw  data. 
Where  missing  data  were  found,    linearly   interpolated  values  were  substituted. 
Since  the  points   in  Figure   7,   unlike  those  in  Figure   3,   are  not   connected,   the 
quantization  of  the  data   (0.1C)  can  be  clearly  seen. 

The  most   interesting  feature  of  Figure  7  is   a  definite  upward  average 
tilt  to  the  data.     A   least-squares   fit  of  the  linear  trend  equation 
y(t)   =  a  +  6t     to  the  data   gives   the  fit     y(t)=10.9  +  0.00037t,   where     t   is   in 
days,   and  y(t)   is  water  temperature  in  °C     on  day  t.     This   indicates   an 
"average"  temperature  rising  from  10.9  to  12.52C  over  the  twelve  year  measure- 
ment period   (t   =  4380).    Note  that  this    linear  trend   is   simply   a  convenient  way 
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of   representing  the  overall    increase  in  temperature  over  the  12-year  period. 
Higher  order  polynomial   fits  provide  only   slight   improvement.      We  may,    how- 
ever,  be  observing  only  a  part   of  a  much   longer  cycle  of  fixed  or   variable 
period  which  would  be  revealed  in  a  much   longer  data  series.   The  statistical 
significance  of  the  trend  parameter       6  =0.00037     (i.e.   whether  it   is   really 
different  from  zero)  is  complicated  by  the  presence  of  cycles  and  possible 
short  term  correlation  in  the  data.     Thus  standard  tests  that  6=0  based  on 
independent  data  are  invalid.   However,   using  simulation  techniques  described 
later,  we  have  found  the  standard  deviations  of  6  to  be  0.000011  and  hence 
reject  the  hypothesis  that  6  =  0  at  the  95%  level   of  confidence. 

The  gross  features  of  the  Granite  Canyon  SST's   clearly  include  annual 
cycles  and  it   is   apparent  from  Figures  2  and  7  that  these  annual   cycles 
contain  inter-annual    variability  as  well.    In  particular  the  oceanic     El   Ninos 
of  1972,   1976,   and  1983  are  clearly  visable.     The  possibility  of  a  temperature 
anomaly   in  1979  is  also  visable.     Although  the  effects  of  El   Nino     events  are 
often  restricted  to  equatorial    regions  and  to  regions   along  the  South  American 
coast,   their  influence  does  upon  occasion  extend  at   least  as  far  north  as 
central  California   (Bretschneider  and  McLain,    1983). 

Underlying  patterns   in  the  data  can  be  exposed  by  performing  local 
smoothing.      (Since  the  computations   in  this  paper  were  done  using 
the  APL  language   (see  e.g.   Gilman  &  Rose,    1976)  with  an  experimental  APL 
Graphics  package,  GRAFSTAT3,   the  smoothing  was  performed  using  the  program 
FILTER   from  Anscombe   (1981,   p392).)     The  smoothing  computes   a  cosine-weighted 
moving  average  of  extent  m  to  give  a  smoothed  time  series       y(t)    : 
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m-1 
y(t)  =   I       y(t)w(t-j)  ,  (3.1) 

J-(m-l) 

where  the  w(k)«  ^-  (l-cos[2ir(k+l)/(m+l)]),  for  k=0,±l,±2,. .  ,±(m-l),   (3.2) 


are  equal -spaced  ordinates  of  l-cos(2t).  The  gain  of  this  filter  is  given  by 
r,    %       sin  [w(m+l)2] 


u)(m+l)[l-a>2(m+l)2(2Tr)2]  (3*3) 

Details   are  given  in  Anscombe   (1981,   p. 156). 

Smoothings  of  the  Granite  Canyon  SST's  for  various  values  of   (2m+l) 
together  with  the  number  of  filter  weights  used  are  given  in  Figure  8.      It   is 
clear  that  61  weights    (or  81  weights)   resolve  the  gross  details  adequately. 
It  is   interesting  to  note  that  there  is   not  only  a  sharp  drop   in  temperature 
in  the  spring,    but  also  a  rather  clearly  defined  and  relatively  abrupt  rise  in 
temperature  in  the  fall.   We  note  that  both  the  spring  and  fall   transitions 
occur  over  periods   relatively  short   in  comparison  with  the  annual   cycle. 

Since  well-defined  spring  transitions  could  be  clearly  identified  in  at 
least  6  out  of  the  12  years,    it  was  decided  to  tabulate  these  events  with 
respect  to  their  dates  of  occurrence,   duration,   and  associated  changes   in 
temperature.     The  results  are  shown  in  Table  1  below.     Because  the  spring 
transitions  add  complexity  to  the  annual   temperature  cycle,   the  harmonic 
content   of  this   data  is  undoubtedly   richer  by   their  presence.      We  will    say 
more  about  the  annual   variation  in  temperature  at  Granite  Canyon   in  Section 
10. 
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Table  1* 
YEAR   STARTING  DATE    DURATION     DELTA  T       dT/dt 


1973 

March  9 

7  days 

13.0  ■ 

-  10. OC 

0.43°/Day 

1974 

March  31 

5  days 

12.3  ■ 

-  10.0 

0.46 

1977 

Feb.  22 

10  days 

13.1  ■ 

-  9.3 

0.38 

1978 

April  25 

9  days 

13.0  ■ 

-  8.5 

0.50 

1980 

March  11 

7  days 

14.5  ■ 

-  10.4 

0.59 

1981 

March  23 

8  days 

13.1  ■ 

-  10.0 

0.39 

*  A  lesser  transition  was  observed  on  December  19,  of  1981.  Whether  or  not 
this  event  corresponds  to  a  "spring"  transition  is  not  known. 


A  detailed  statistical  analysis  of  the  time  series  of  SSTs  will  now  be 
given.  This  will  be  used  in  Sections  8  and  9  to  give  a  more  detailed  exam- 
ination, statistical  and  oceanographic,  of  the  SSTs. 

4.   Cyclic  fit  and  decomposition  of  the  data:  General  ideas 

We  proceed  to  decompose  the  SST  time  series  into  an  additive  model  of  the 

form 

y(t)  =  m(t)  +  s(t)  +  e(t),  (4.1) 

where  the  trend  m(t)  can  be  defined  loosely  as  a  long  term  change  in  the  mean 
(Chatfield,  1980,  p. 13),  s(t)  consists  of  seasonal  and  cyclic  changes,  and 
e(t)  is  a  stationary  random  sequence  which  describes  the  irregular  fluctua- 
tions. The  e(t)  sequence  is  assumed  to  have  mean  zero,  (i.e.  E(e(t)J=  0), 
constant  variance,  and  possibly  some  serial  dependence.  That  m(t),  s(t)  and 
e(t)  are  treated  as  additive  implies  the  gross  assumption  that  this  dependency 
structure  in  the  e(t),  as  well  as  the  variance,  is  not  affected  by  m(t),  or 
s(t). 
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The  model  (4.1)  is  essentially  a  conceptual  tool  for  thinking  with  and 
its  usefulness  and  approximate  validity  must  be  established  from  the  data. 
Thus  to  define  "scales"  in  the  data  these  different  components  must  be 
identified. 

The  term  m(t)  was  discussed  in  Section  3  and  here  we  let 

y:(t)  =  y(t)  -  m(t)  =  s(t)  +  e(t)  (4.2) 

insofar  as   it   can  be  observed.      With   the  estimate     m(t)   =   10.9  +  0.00U374t 
removed  from    y(t)     we  denote  the  resultant   as  the  first   residual,    r\(t), 
ignoring  for  the  present  possible  interaction   between  the  estimation  of  m(t), 
s(t)  and  e(t). 

In  this  section  we  estimate  the  structure  of  s(t)  and  e(t)  by  analyzing 
r\  (t )   .     Note  that     s(t)     includes  definite  seasonal   terms,   such 
as   the  obvious   annual   effect   in     E(y(t))    (and   its   harmonics),    and  terms  with 
periods   longer  than  a  year  which  might  be  interpreted  as   global   fits  to  the 
obvious   long-term  modulation  of  the  annual   cycle  which   can   be  seen   in  figures 
2  and  6.         Again  there  may   be  cycles   at  shorter  periods   associated  with   the 
tides  and  winds;   these  are  generally  difficult  to  distinguish  from 
quasi-cycles   "generated"  by  highly  correlated  errors  e(t).     Thus  we  assume 
that 


J 

s(t)   =     I     (y.-  cosU.t)  +  B-   sin(a).t)}  (4.3) 

,■  _  j       J  J  J  J 


J 

I       A.   cosU.t   +  6   )    ,  (4.4) 

•j  =  i       J  J  J 
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where     w.=  2irf  .   =  2tt/T.   are   (unknown)   frequencies   corresponding  to  periods     T., 
j  j  j  j 

221/2  -1 

A.   =   (y.   +  3.)  is     the  amplitude,        8.=  tan     (-6/y-)     is   the  phase,   and 

J     is   the  number  of  cycles  which  comprise     s(t).     Again  we  note  that   cycles 

with  periods   longer  than   1  year  may   represent   long-term  random  fluctuations. 

The  representation    (4.3)  of  the  component  s(t)  of  E(y(t))  as  a  sum  of 

deterministic  cycles  with  fixed  frequencies   is  a  convenience  only.      In 

particular  it  might   be  hazardous  to  predict   very   long  term  behavior  on  the 

basis  of  this  model   for  s(t)  and  m(t).     Also,  while  our  model   may  identify 

cycles  with  periods  shorter  than  one  year,   such  cycles  may  not   represent   real 

cyclic  behavior,   but   rather  they  may  serve  only  to  accommodate  the 

non-sinusoidal   nature  of  the  year  cycle  and  the  very  sharp  transition  in  that 

cycle  that  often  occurs  during  spring   (see  Table  1).   Because  our  data  extends 

for  only   12  years,    it  will    be   impossible  to  detect  any  cycles  with  periods 

longer  than  half  that   length,   or  6  years. 

To  find   "hidden  periodicities"  in  s(t)   is   a  formidable  task,   even  if  the 

spectrum  of  the  random  component  e(t),    namely 

00 

f (co)   =  —  {1  +  2     I     p(k)   cos(kai)}    ,  0  <  a)  _<  it   ,  (4.5) 

11  k=l 

where  p(j)   =  corr(e(t+j )e(t )),  were  constant    (flat)   i.e.   the  e(t)  are   "white 
noise".    (This   is  the  classical   problem  of   "hidden  periodicities"  going  back  to 
the  original   work   of  Schuster   (1898)   and  Fisher   (1929).)     To  find  these 
periodicities   in  the  presence  of  a  non-flat  spectrum  for  the  e(t)'s,  we  proceed 
as   follows: 

(1)  Examine  the  periodogram  of  the  first  residuals,  r,(t),  to  find 

possibly  significant  periodicities.     Estimate     s(t)  as     s(t). 

Given  the  chosen  frequencies  to.,  this  involves  obtaining  least 
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squares   estimates   of  the     y    '-s   and  B-'s      in  the   linear  model    (4.3). 
The  least  squares   estimates   of  the  amplitudes     A.     are  directly  pro- 
portional  to  the  periodogram  values   at  the  given  frequencies. 

(2)     Examine  the  second  residuals     r2(t )=y (t )-s (t )-m(t )=r^(t )-s (t ) 

to  determine  the  properties   of     e(t).      If     e(t)     can   re   represented, 
say,   as  a  pth-order  autoregressi ve  process, 

e(t)    =  a1e(t-l)   +  ...    +  a   e(t-p)   +  n(t),  (4.6) 

where  the     n(t)   are  uncorrected.    Then     we  estimate     a,      through     a 
by  least  squares  to  be     a,     through     a       and     "pre-whiten"  the  first 
residuals,     r,(t),   as 

r3(t)  =  rx(t)   -  a^^t-1)   ...    -a   r^t-p).  (4.7) 

We  then  examine  the  periodogram  of  these  third  residuals   to 
ascertain  whether  the     s(t)     chosen  in  step    (1)   is  still   a 
"reasonable"  representation  of  the  seasonal   and  cyclic  changes 
(global    fit)   for  the  data.   The   "prewhitening"   is  performed  because 
highly  autocorrelated  data  will   elevate  spectral    values   near  zero 
frequency,   making  it  difficult  to  decide  whether  a  given  spectral 
peak   is  due  to  the  correlation  or  due  to  a  cycle  in  the  data.     Once 
the  autocorrelated  component   is   removed,  we  obtain  a  clearer  picture 
of  what  cycles  should  actually     be  included  in     s(t). 

(3)     Using     To(t)     to  point  out  a  possibly   refined     s(t),      re-examine  the 
structure  of     e(t)   via  the   (possibly   refined)  second  residuals, 

r~(t  )=y  (t  )-m(t)-s(t ).      Knowledge  of  the  structure  of     e(t)   is 
important   for  determing   "local   scales"  and  obtaining  the   "whitened" 
residuals     r«(t).     Thus,    if   it   is  determined  that     e(t)  has     pth 
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order  autoregressi ve  structure,   possibly  different  from  that   in   (2) 
above,   compute  fourth   residuals 

r4(t)   =  r2(t)   -  a^t-l)    ...    -a   r2(t-p).  (4.8) 

These  residuals  are  used  firstly  to  establish  the  overall  goodness 
of  fit  of  the  estimated  model,  since  the  r*(t)  process  will  have 
a  flat  spectrum  and  a  zero  level  correlogram,  (i.e.,  "white  noise") 
if  the  model  is  reasonable.  Secondly  they  can  be  used  (e.g.  Tsiao  and 
Box,  1982)  to  establish  cross-correlations  of  the  Granite  Canyon  data 
with  other  SST  data.  Thirdly,  they  provide  an  estimated  distribution 
of  the  n(t)  which  will  be  useful  later  for  modelling  and  for  checking 
assumptions  during  statistical  testing. 

5.   Global  fit  and  cyclic  components 

Figure  11  shows,  in  the  top  panel,  a  cyclic  fit  s(t)  with  J  =  13  terms 
to  the  first  residuals  (original  data  minus  long  term  linear  trend  m(t)).  The 
data,  since  leap  days  have  been  removed,  consists  of  N  =  4380  daily  observations 
of  temperature.  Thus  if  angular  frequency  is  re-expressed  as  discrete  frequencies 
of  the  form  w  =  2-iTp/N,  where  p  is  the  "number  of  cycles  per  period  of  4381) 
days",  then  the  cycles  removed  from  the  data  as  components  of  s(t),  have 
frequencies  p  =  12,24,94,  corresponding  to  periods  of  1  year,  6  months  and 
approximately  1.5  months,  with  the  remaining  p  =  2,3,4,5,6,7,8,9,11,18. 

The  determination  of  these  frequencies  follows  the  steps  in  Section  4. 
First,  Figure  10  shows  the  estimated  serial  correlations  of  the  first 

residuals 

N-k 

P   (k)  =  I  (y(t)  -m(t))  (y(t+k)-m(t))/{(N-k)S2}       (5.1) 
t=l 

N-k 
=  I     r1(t)r,(t+k)   {(N-k)S2} 
t=l  L         l 
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where  S2  is  the  sample  variance  of  the  r  (t)'s.  This  correlogram  is 

dominated  by  the  yearly  cycle  and  therefore  can  not  be  used  to  postulate  a 

model  for  the  e(t),  or  to  identify  more  than  the  yearly  cycle  itself. 

A  harmonic  spectral  analysis  of  the  p(k)'s  is  given  by  the  periodogram 

of  the  first  residuals,  defined  as 

In(u  )  =  {A*  (w  )  +  B*  (o)p)}  /(2irN),   p  =  l,..,N/2  ,       (5.2) 

where 

N  N 

AN(o)p)=  I     rx(t)  cos(tcop),  and  BN(wp)  =  I     rx(t)  sin(tu>p), 

for  p  =  l,...,N/2,  are  the  real  and  complex  parts  of  the  discrete  Fourier 
transform  of  r,(t),  t=l,  ..,  N.  Here  w  =  -rp-  .  Since  the  expected  value  of 

2 
the  periodogram  is  (asymptotically),  for  uncorrelated  data,  a   /2tt,  we  use 

a  normalized  periodogram.  Thus  the  normalized  periodogram  is  defined  to  be 
the  result  (5.2)  multiplied  by  2-n/S2,  where  S   is  the  sample  variance  of  the 
data  being  considered,  namely  r\(t). 

The  logarithm  of  the  normalized  periodogram  of  the  first  residuals  from 
the  Granite  Canyon  data  is  shown  in  Figure  9  with  confidence  bands  for 
individual  components  and  for  the  maximum  value  of  the  periodogram  components. 
These  are  based  on  the  approximation  that  each  value  in  the  normalized  period- 
ogram is  independently  distributed  as  a  mean  one  exponential  random  variable 
(chi -square  with  2  degrees  of  freedom).  The  dominant  frequency  is  p  =  12, 
i.e.,  a  period  of  one  year.  There  is,  in  the  expanded  scale  periodogram  in  the 
top  panel  of  Figure  9,  a  significant  harmonic  at  p  =  24  and  at  p  =  94.  This 
latter  frequency  is  surprising,  since  it  is  not  exactly  a  harmonic  of  p  =  12. 
The  full  periodogram  in  the  lower  panel  suggests,  even  though  it  is  unsmoothed, 
an  autoregressi ve  type  spectrum  for  the  e(t).  Thus,  determination  of  cycles  at 
frequencies  below  12  is  difficult.  Proceeding  to  the  prewhitening,  frequencies 
p  =  2,3,4,5,6,7,8,9,11  and  18  are  selected  for  removal,  in  addition  to  p  =  12, 
24  and  94. 
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Removal  of  these  frequencies  gives  the  second  residuals  whose  periodogram 
is  shown  in  Figure  12.  On  the  expanded  scale,  note  the  zeros  where  cycles  have 
been  removed;  this  is  because  the  least  squares  estimates  of  the  amplitudes  ir 
(4.4)  are  exactly  proportional  to  the  corresponding  periodogram  values.  Again 
note  in  the  lower  panel  of  Figure  12  the  autoregressi ve-1 ike  spectrum.  (A 
spectrum  for  a  theoretical  AR(2)  process  using  parameters  computed  below  is 
shown  for  comparison.) 

Figure  13  shows  the  estimated  autocorrelations  of  the  second  residuals  as 
well  as  the  first  5U  correlations  for  four  theoretical  AR(P)  models  computed 
from  the  first  four  estimated  correlations,  via  the  Yule-Walker  equations 
(Chatfield,  1980,  p. 48).  Note  in  Figure  13  that  the  AR(P)models  of  order 
greater  than  2  do  not  give  a  much  better  fit  than  the  AR(2)  model  to  the 
estimated  autocorrelations. 

Notice  also  that  the  persistent  hump  at  about  21  days.  As  yet  we  have  no 
explanation  for  this  effect.  The  autocorrelations  in  Figure  12  are  extended  to 
show  lags  beyond  50  in  Figure  13,  but  there  is  no  indication  of  a  cycle  of 
period  21.  Such  a  cycle  would  show  up  as  a  wave  in  the  correlogram  with  period 
21,  not  simply  as  a  peak  at  lag  21  days. 


Table  3 
Estimated  correlation  coefficients  for  the  residuals  r  (t)  and  estimatec 
coefficients  for  autoregressi ve  models  of  orders  p  =  1,2,3  and  4.  Note  that 
the  underlined  values  are  the  first  four  partial  correlations  for  the  data. 


Correlations 


P(D 
5(2) 
P(3) 

P(4) 


0.8511 
0.6973 
0.5718 
0.4710 


Estimated  coefficients 


ax   =  0.8511 
a:   =  0.9343;    a. 


■0.0975 


al   =  0.9355;   a2   =  0.1098;   a3   =  0.0129 
a:   =  0.955;     a2  =  0.109;     a3  =  0.007;   al 


=   0.0006 
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Consequently  the  first   residuals   are  pre-whitened  as   in   (4.7)   above,   using 
p   =  2  and     a\    =  0.9343  and     a2   =  -0.0975  derived  from  the  values    *p(l)   =  0.8511 
and     p(2)   =  0.6973.     The  resulting  third  residuals,      r3(t),    and  their  auto- 
correlations  are  shown  in  Figure  15  with  the  logarithm  of  the  normalized 
periodogram  in  Figure   16.     The  choice  of  frequencies   removed  on  the  basis   of 
the  first  residuals   is  confirmed.     To  explain  the  hump   in  the  autocorrelogram 
(Figure   14),   several   other  cyclic  components  were   removed,   e.g.    those  at 
p   =  19,   20  and  243.     These  changed  neither  the  autocorrelations   nor  the 
residual    variance  of  the  data. 

We  discuss  now  details  of  the  detrending  of  the  data.     Table  4  shows  the 
sample  variance  of  the  original   data   (raw  data  minus  grand  mean),   the  first 
residuals    (removal   of  linear  trend),   the  second  stage  (removal   of  linear  trend 
and     s(t)     to  give     r   (t)),   and  the  fourth   residuals  which  will    be  discussed 
below. 


Tabl 

e  4 

Stage 

Samp 

le  Variance 

Sa 

mple  s.d. 

%  variance   removed 

Raw  data 

2.624 

1.62 

8.4% 

rx(t) 

2.403 

1.55 

51.9% 

r2(t) 

1.040 

1.02 

28.7% 

r4(t) 

0.288 

0.537 

10.9%     Remainder 

Note  that  about  half  of  the  variance  is  accounted  for  by  the  cyclic 
component  s(t).     A  break-down  for  the  J   =   13  cyclic  components   is   given  below 
in  Table  5.     Components  are  listed  in  descending  order  of  magnitude  of  their 
amplitudes,     A.   .     The  variance  after  removal   of  each  of  these   (orthogonal) 

u 

components   is   given  in  the  last  column.     The  annual   cycle  accounts  for  a  large 
part  of  the  variation,   namely   (. 81)/1. 36  x  100%  =  59.6%     of  this  variance. 

These  estimated  coefficients  will   be  used  in  later  sections   in  a  simula- 
tion study  of  the  model. 

28 


go 
< 

o 

LxJ 

cr 


o 
o 

Ld 
CO 

a: 
o 


o 


UJ 

q: 
o 
o 
O 

< 


1 
1 

f; 

j 

t 

I 

t- 

il 

ij 

I ...I-.... 

S0V1 

i  / 

100 

y 

i 

I 

i 

i 

i 

i 

i 

n 

i 

i 

i 
i 

/                  i- 
1             1             '             I'           J 1              1 

i 

__L 

g 


o 

o 


rO  90  *0  20 

iN3ioiii300  NOiiyisaaoDoinv 


80  90  ►  0  I'O 

1N3OUJ303  NOiHTGdaODOinv 


29 


CM 

< 

i i 

o 


UJ 


I 

Ld 

cl 


O 
u_ 

< 

9 

CD 
UJ 

q: 
q 

X 


1 

z 

2 


in 

< 
Q 


in 


O 

M 


2  C  Z-  «■- 


Jl'O  800  K)0  0 

iN3i3'.dJ303  no:  L\n3aaOD0inv 


30 


to 
< 

9  ^ 

CO  CM 

DC   9? 
<. 

Q     I 

q:  in 

K  3 
L^   9 

o  oo 

-^  Id 
1* 

or  h- 
o  to 
O  £ 

cr 

LjJ 
Q. 


CO 

UJ 

—I 


V 

o 


o 


CO 


o 

> 

u 

VC 

H 

z. 

W 

>- 
u 

g 

z 

u 

UJ 

t-i 

O 

(J- 

O 

LxJ 

a 

2  0  3-  V- 

wvyooaoiy3d  ooi 


►  0  V- 

wvijooaoiy3d  jo  ooi 


31 


Table  5 

j 

PJ 

T.  (years) 

YJ 

Bo 

A. 
J 

6  . 
J 

Variance 
2.40 

1 

12 

1.00 

-0.465 

-1.187 

1.275 

111.4° 

1.59 

2 

6 

2.00 

0.273 

-0.295 

0.402 

47.2° 

1.51 

3 

24 

0.50 

0.342 

-0.162 

0.378 

25.4° 

1.44 

4 

7 

1.71 

0.345 

0.112 

0.363 

342.0° 

1.37 

5 

3 

4.00 

-0.266 

0.211 

0.340 

218.4° 

1.32 

6 

9 

1.33 

-0.101 

0.302 

0.318 

251.6° 

1.28 

7 

2 

6.00 

0.074 

0.275 

0.285 

285.1° 

1.23 

8 

11 

1.09 

0.097 

0.250 

0.268 

291.2° 

1.19 

9 

18 

0.67 

-0.187 

-0.180 

0.260 

136.1° 

1.15 

10 

8 

1.50 

0.079 

-0.243 

0.256 

72.0° 

1.12 

11 

5 

2.40 

-0.176 

-0.174 

0.247 

135.3° 

1.08 

12 

94 

0.13 

-0.155 

-0.152 

0.218 

135.6° 

1.06 

13 

4 

3.00 

0.194 

-0.097 

0.217 

26.6° 

1.04 

6.  Examination 

of  the  random  component,  e(t) 

The  second  residuals,  the  original  data  minus  the  linear  trend  and  the 
seasonal  trend  s(t)  containing  13  components,  are  whitened  (4.8)  under  the 
assumption  of  the  adequacy  of  an  AR(2)  model  to  give  the  fourth  residuals 

r4(t)  =  r:(t)  -  a1r2(t-l)  -  a2r2(t-2). 

The  first  use  of  these  residuals  is  as  a  test  of  the  adequacy  of  the 
modelling  which  has  been  done  up  to  this  point.  There  are  two  issues  which 
arise;  (1)  the  physical  adequacy  of  the  model, and  (2)  the  statistical  adequacy. 
Any  given  model  may  satisfy  either,  neither  or  both  of  these  requirements. 
For  physical  adequacy  we  require  that  the  model  describe  a  substantial  amount 
of  the  variation  in  the  data  with  as  few  terms  as  possible  (i.e.,  the  principle 
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of  parsimony),  and  hope  that  at  least  the  dominant  terms  have  oceanographic 
interpretations.  For  statistical  adequacy  we  have  made  the  assumption  that  the 
fourth  residuals  are  estimates  of  the  n(t)'s  and,  as  such,  should,  if  the  model 
is  reasonable,  be  approximately  uncorrelated. 

The  residuals  themselves  are  plotted  in  the  top  panel  of  Figure  17.  There 
is  an  anomaly  at  about  3100  days  associated  with  isolated  extremes  in  SSTs 
encountered  during  the  fall  of  1979.  It  is  in  fact  a  reminder  that  the  model  we 
have  derived  is  a  global  model  which  will  not  account  for  transients  or,  say, 
subtle  seasonal  effects.  Except  for  that  anomaly,  the  fourth  residuals 
generally  fall  in  the  range  of  ±  2.0  which  is  considerably  smaller  than  the 
range  of  the  original  data  (8.0  to  17.0).  Moreover,  in  Table  4  we  have  shown 
that  the  fourth  residuals,  which  represent  the  unpredictable  part  of  the  data, 
account  for  only  10.1%  of  the  total  variation.  Oceanographic  interpretations 
for  some  of  the  parameters  used  in  constructing  the  model  are  presented  later, 
as  well  as  simulation  results  which  show  how  well  the  model  mimics  in  behavior 
the  data. 

The  autocorrelation  function  for  the  fourth  residuals  is  given  in  the 
lower  panel  of  Figure  17.  The  bands  are  approximate  9b%  confidence  bands  for 
individual  correlations.  Again  there  may  be  an  anomaly  at  lag  20,  but  since 
this  is  a  maximum  of  100  estimated  correlations  it  probably  is  not 
significantly  large.  A  more  precise  test  that  the  r  (t)s  are  uncorrelated 
is  given  in  Figure  18. 

The  top  panel  is  the  unsmoothed,  normalized  periodogram  and  the 
bottom  panel  is  the  cumulative  periodogram.  Barely  distinguishable  upper 
and  lower  95%  and  99%  confidence  intervals  for  a  flat  spectrum  based  on 
Kolmogorov-Smi rnov  (KS)  statistics  are  shown.  Since  the  cumulative  periodogram 
lies  within  these  bands,  the  Durbin  test  for  a  flat  spectrum  (see  Jenkins  and 
Watts,  1969,  p. 234-7,  or  Cox  and  Lewis,  1966,  p. 169)  accepts  the  hypothesis  of 
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a  flat  spectrum.     The  value  of  the  KS  statistic  is   given  as  0.5136  and  compares 
to  the  upper  5%  point  of  the  statistic  which  has   value  1.355.    In  addition  the 
Anderson-Darling  test  statistic   (value  0.9139)   is   given;   this  test  statistic  is 
more  sensitive  to  departures   at   low  or  high  frequencies    (Cox  and  Lewis,   1966, 
p. 169)  and  in  this  case  accepts  the  hypothesis  of  no  correlation. 

Thus,  we  have  established  the  adequacy  of  the  model.     Note  that  if  the 
e(t)  process   is  a  Normal   process,   then  a  flat  spectrum  implies   independence, 
but  not  otherwise.     For  this   reason,   and  in  order  to  complete  the  model,  we 
examine  the  distribution  of  the  fourth  residuals. 

Figure  19  shows,   in  the  top  panel,   a  crude  histogram  of  the  fourth 
residuals,   together  with  sample  statistics.     The  estimated  kurtosis  value  of 
3.61,  which  for  the  sample  size  of     N  =  4380     is  significantly  different  from 
zero,    indicates  non-normality   in  the  residuals.     This   is  confirmed  in  the 
probability  plot   in  the  lower  panel.     This   is  a  plot  of  the  ordered 
residuals,   say     X,.s,   against  their  approximate  expected  value     $     (i/N+1) 

where     $     (.)     is  the  inverse  normal   probability  function.     The  curvature  in 
this  plot  and  the  outliers   indicate  that  the  distribution  of  the     e(t)'s   is 
somewhat  non-normal,    in  fact  pinched  in  the  middle  and  long  in  the  tails.     The 
distribution  appears  symmetric;  this   is   reflected  in  the  estimated  coefficient 
of  skewness  of     -0.00147. 

No  attempt  has  been  made  to  model   these     e(t)'s   in  more  detail.     Their 
standard  deviation  is,   from  the  last  line  of  Table  4,   estimated  to  be  0.537. 
The  autoregressi ve  parameters,     a     and  a   ,     for  the  AR(2)  process  are  estimated 
in  Table  3  to  be  0.9343  and  -0.0975  respectively  and  the  coefficients  of     s(t) 
are  given  in  Table  5.     The  long-term  trend  is     m(t)   =  10.9  +  0.000374t. 
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7.   Comments  on  the  model  and  the  modelling 

In  Figure  8  we  have  given  a  local  smoothing  of  the  data  which  shows  the 
main,  gross  features  of  the  SSTs.  Although  there  are  many  ways  in  which  the 
data  can  be  smoothed  and  there  is  the  choice  of  an  appropriate  bandwidth,  the 
procedure  is  free  of  the  need  to  postulate  a  global  model.  In  fact  the 
smoothing  (Figure  8)  and  the  global  fit  of  s(t)  in  Figure  11  give  very 
similar  information.  Why  then  consider  global  models  of  m(t)  and  s(t)  and  the 
subsequent  modelling  of  e(t)? 

The  answer  lies  in  our  objectives. 

( 1 )  We  want  to  use  the  model  to  identify  relevant  "scales"  in  the  data? 
The  question  of  scales  is  rather  nebulous,  with  different  answers 
coming  out  of  different  parts  of  the  model.  With  regard  to  s(t), 
there  is  clearly  a  "scale"  of  one  year  and  longer  scales  that  apply, 
for  example,  to  the  occurrences  of  El  Ninos.  Locally,  the  high 
correlation  in  e(t)  gives  another  scale,  measured  perhaps  by  the 
autocorrelation  function  of  the  e(t)'s  given  in  Figure  14.  And,  of 
course,  if  diurnal  effects  are  of  interest,  then  there  is  a  scale  of 
interest  shorter  than  the  one  day  sampling  interval  used  to  collect 
the  data. 

(2)  With  respect  to  scales,  we  want  to  know  how  often  should  the  data  be 
sampled?  This  is  a  question  which  cannot  be  answered  without  a 
specific  objective  in  mind.  For  instance,  if  we  want  information 
about  daily  variation,  we  need  to  sample  more  often  than  once  a  day. 
If  we  are  concerned  with  longer-term  variation,  we  can  relate  the 
sampling  rate  to  the  accuracy  of  predictions.  One  component  of  the 
prediction  is  the  cyclic  and  long-term  trend;  how  accurately  are  they 
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estimated  and  what  would  be  lost  by  sampling  once  every  two  or  three 
days?  In  Section  9  we  examine  sampling  variation  in  the  parameters 
in  m(t)  and  s(t)  by  generating  (simulating)  independent  realization 
of  SST  based  on  the  model. 

(3)  The  model  can  be  used  to  interpolate  missing  data.  This  is  not  a 
major  problem  in  the  Granite  Canyon  data  set  and  thus  it  will  not  be 
discussed.  However,  gaps  are  a  frequent  problem  with  many 
oceanographic  data  sets  (e.g.  in  the  Farallon  Island  data  mentioned 
previously). 

( 4 )  The  model  can  be  used  to  predi ct  or  forecast  future  temperatures, 
i.e.  given  the  data  up  to  today,  what  will  they  be  for  several 
subsequent  days.  We  show  in  the  next  section  how  the  model  is  used 
to  give  predictions  for  up  to  two  days  in  advance  of  current  data. 
These  predictions  will  be  shown  to  be  better  than  the  usual 
climatological  prediction  based  on  averages. 

If  we  attempt  to  use  a  climatological  model  (i.e.  an  average 
across  the  12  years  or  an  average  across  the  smoothed  values  for  12 
years)  to  address  the  questions  above,  we  sacrifice: 

(1)  The  variation  from  year-to-year,  so  that  events  like 

El  Ninos  and  the  long  term  trend  represented  by  m(t)  are 
masked  by  the  averaging. 

(2)  The  local  information  represented  by  e(t).   It 
seems  intuitively  reasonable  to  make  use  of  the  water 
temperature  today  if  we  need  to  make  a  prediction  for 
tomorrow.  Quantitatively,  Table  4  has  shown  that  this 
local  information  accounts  for  28.7%  of  the  variation  in  the 
data. 
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(3)  Decomposition  of  the  different  factors  affecting  sea-surface 
temperature,  retaining  only  a  picture  of  the  gross  behavior. 
Also  we  can  no  longer  consider  separately  the  yearly,  6 
monthly,  or  tidal  cycles. 
Hence,  we  find  the  model  to  be  more  useful  than  climatology  for 
summarizing  the  data.  We  use  the  next  two  sections  to  demonstrate  its  value. 

8.   Forecast  of  SSTs  from  March  1,  1983  to  March  30,  1983 

We  now  forecast  SSTs,  one  and  two  days  ahead,  using  the  model  developed  in 
the  previous  sections.  Since  data  are  now  available  for  the  month  of  March, 
1983,  a  period  beyond  the  12  year  period  used  in  constructing  the  model, 
predictions  are  made  for  comparison  with  the  observations.  The  previous  twelve 
years  of  data  are  used  to  estimate  the  parameters  in  the  model.  These  forecast 
values  are  compared  to  the  actual  observations  and  to  a  climatological  forecast 
using  the  average  temperature  for  each  given  day  across  the  twelve  year 
ensemble. 

Given  an  autoregressi ve  process,  i.e,  the  second  residuals  r  (t)  as  in 
equation  (4.1)  of  Section  4,  with  known  coefficients  a  and  a  , 

r2(t)  =  a1r2(t-l)  +  a2r2(t-2)  +  e(t), 

and  that  the  i.i.d  residuals  are  normally  distributed  with  mean  zero,  then  the 
best  predictors  one  and  two  steps  ahead  are   obtained  by  setting  e(t+2)  and 
e(t+l)  equal  to  their  expected  values  of  zero.  We  obtain 

Pred  (r2(t+l))  =  a^t)  +  a2r2(t-l)  , 
and  2  * 


Pred  (r2(t+2))  =  (a^  +a2)r2(t )  +  a^r^t-l)  , 
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where     a     and  a       are   least  squares   estimates   of     a     and     a     based  on   the 
observed   12  years   of  data.      For  periods   far  enough   ahead,    these  predicted 
values  would  eventually  converge  to  zero.      This    is   simpler  to  see  if     a     =  0 
so  that  we  have  an  AR(1)   process.     Then     Pred(r2(t+k )  J   =  a-,    i~2(t)   . 

We  also  have  estimates  rfi(t)  and  s(t)  for  the  linear  and  seasonal  trends 
for  t=l  to  4380,  and  these  can  be  computed  beyond  that  range  to  give 
projections   m(t+l),  m(t+2),  and  s(t+l),  s(t+2).  The  forecasts  are  then 

Pred(y(t+1))  =  m(t  +  l)  +  s(t+l)  +  Pred(r2(t+1) } 

Pred(y(t+2))  =  m(t+2)  +  s(t+2)  +  Pred(r2(t+2) ) 

These  forecasts  have  been  computed  for  the  moving  twelve  year  period  out 
to  March  30,  1983,  together  with  the  cl imatological  forecast,  and  all  are  given 
in  Table  6. 
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Date 


Table  6 

Forecasts   of  one  month's  SSTs   at  Granite  Canyon 

One-step 
One-step  Two-step  Climatol. 

Actual  y(t)   Forecast     Error     Forecast  Error  Forecast  Error 


March 

1. 

1983 

4381 

14.3 

14.0 

-0.3 

11.5 

-2.8 

4382 

14.4 

14.0 

-0.4 

13.7 

-0.7 

11.4 

-3.0 

4383 

14.2 

14.1 

-0.1 

13.7 

-0.5 

11.5 

-2.7 

4384 

14.5 

13.9 

+0.6 

13.8 

-0.7 

11.4 

-3.1 

4385 

14.3 

14.2 

-0.1 

13.6 

-0.7 

11.2 

-3.1 

4386 

14.4 

14.0 

-0.4 

13.8 

-0.6 

11.3 

-3.1 

4387 

14.2 

14.1 

-0.1 

13.7 

-0.5 

11.4 

-2.8 

March 

8, 

1983 

4388 

14.8 

13.9 

-0.9 

13.8 

-1.0 

11.6 

-3.2 

4389 

14.9 

14.5 

-0.4 

13.6 

-1.3 

11.6 

-3.3 

4390 

14.7 

14.5 

-0.2 

14.1 

-0.8 

11.5 

-3.2 

4391 

14.8 

14.3 

-0.5 

14.1 

-0.7 

11.4 

-3.4 

4392 

14.0 

14.4 

+0.4 

14.0 

0.0 

11.6 

-2.4 

4393 

14.2 

13.7 

-0.5 

14.1 

-0.1 

11.1 

-3.1 

March 

14, 

1983 

4394 

14.5 

13.9 

-0.6 

13.5 

-1.0 

11.0 

-3.5 

4395 

13.6 

14.2 

+0.6 

13.7 

+0.1 

10.7 

-2.9 

4396 

13.5 

13.3 

-0.2 

13.9 

-0.4 

11.0 

-2.5 

4397 

13.8 

13.3 

-0.5 

13.2 

-0.6 

10.9 

-1.9 

4398 

13.7 

13.6 

-0.1 

13.2 

-0.5 

10.9 

-1.8 

4399 

13.7 

13.5 

-0.2 

13.4 

-0.3 

11.0 

-2.7 

4400 

13.7 

13.5 

-0.2 

13.3 

-0.4 

10.9 

-2.8 

March 

21, 

1983 

4401 

13.5 

13.5 

0.0 

13.4 

-0.1 

11.1 

-2.4 

4402 

13.8 

13.3 

-0.5 

13.4 

-0.4 

11.0 

-2.8 

4403 

13.2 

13.6 

+0.4 

13.2 

0.0 

11.0 

-2.2 

4404 

13.5 

13.0 

+0.5 

13.4 

-0.1 

11.0 

-2.5 

4405 

13.3 

13.4 

+0.1 

12.9 

+0.4 

10.8 

-2.5 

4406 

12.7 

13.2 

+0.5 

13.2 

+0.5 

10.6 

-2.1 

4407 

12.6 

12.6 

0.0 

13.0 

+0.4 

10.2 

-2.4 

March 

28, 

1983 

4408 

12.5 

12.6 

+0.1 

12.6 

+0.1 

10.2 

-2.3 

4409 

12.4 

12.5 

+0.1 

12.5 

+0.1 

10.2 

-2.2 

4410 

12.8 

12.4 

-0.4 

12.5 
12.4 

-0.3 

10.5 

-2.3 

Av 

-0.14 

Av 

-.38 

Av 

-2.7 

s.d 

0.37 

s.d 

0.42 

s.d. 

0.45 

m.  s.e 

0.40 

m.  s.  e, 

.  0.57 

m.  s.  e. 

2.74 

First  one  can  see  that  the  one-day-ahead  forecast     Pred(y(t+1))  and  the 
two-day-ahead  forecast,     Pred(y (t+2) ),   are  usually  too  low   (i.e.    negatively 
biased),   the  average  differences  being  -0.14  and  -0.38  respectively.     The 
climatological   forecast   is  even  worse.     Of  course,   this   is  probably  due  to 
the  effect   of  the  current  El   Nino.      The  global   prediction  of  m(t)+s(t)   cannot 
quite  track   the  anomalous  warming,    although  the  autoregressi ve   residuals 
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compensate  to  some  extent.  The  cl imatological  forecast  is  quite  far  off  and  the 
comparison  because  of  the  El  Nino  is  probably  not  fair,  but,  nevertheless,  the 
model  given  in  this  paper  does  give  a  better  forecast.  Note  that  the  s.d.  of  the 
of  the  global  forecast  increases  as  one  goes  from  one  step  ahead  to  two  steps 
ahead,  as  does  the  bias.  This  is  to  be  expected  since  less  local  information  is 
available  for  the  two-day-ahead  prediction.  A  better  over-all  measure  of  the 
goodness  of  the  forecast  is  the  mean  square  error.  This  is  given  in  all  three 
cases  in  the  last  row  of  Table  6. 

9.   Simulation  results 

Here  we  use  the  model  for  the  Granite  Canyon  SSTs  to  generate  independent 
replications  of  the  time  series.  These  can  be  considered  as  members  of  an 
ensemble  from  such  a  process,  with  the  proviso  that  the  model  and  the  model 
parameters  used  to  generate  the  time  series  are   those  estimated  from  the 
original  Granite  Canyon  series. 

One  such  realization  is  shown  in  Figure  20.  The  top  panel  is  the  original 
series  and  the  bottom  panel  is  the  simulated  series.  Note  that  the  simulated 
series  is  generated  out  one  year  beyond  March  1,  1983  and  this  projection  is 
shown  in  Figure  21.  This  is  not  a  true  forecast,  except  insofar  as  the  trend 
m(t)  and  the  cycles  s(t)  account  for  a  major  part  of  the  variation  in  the 
Granite  Canyon  SSTs.  We  note  however,  that  the  one  year  simulation  (1983-84) 
does  contain  an  apparent  spring  transition  and  thus  introduces  a  note  of  reality 
into  the  results. 

Ten  independent  realizations  of  the  Granite  Canyon  SSTs  were  generated  and 
used  to  compute  estimated  parameter  values.   In  these  simulations,  the  "white 
noise"  sequence  n(t)  was  taken  to  be  i.i.d  Normal  random  variables,  despite  the 
departures  from  Normality  shown  in  Figure  19.  An  alternative  would  have  been  to 
"bootstrap"  the  process  by  using  a  random  permutation  of  the  r  (t)'s. 
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In  Table   7  we  give  the  results   of  the  simulation  for  the  parameters   a  and  6 
in  the   long  term  trend,    for  the  serial   correlations   at   lags  one  and  two  in  the 
second  residuals     r   (t),   for  the  associated  coefficients     a     and  a       in  the  AR(2) 
model,    and  for  the  amplitudes     A     and     A2     for  the  two  dominant   frequencies   in 
s(t)  with  periods  of  one  year  and  two  years.     These  are  also  denoted     A     and  A 
in  Table  5.   The  nominal   values  are  the  values  estimated  from  the  real   data  which 
were  used  to  generate  the  simulations. 

Table   7 
Simulation   results   -  sampli  ng  every  poi  nt 
Parameter       Nominal    value  Min  Max  Average s.d. 


a 

10.9 

10.78 

10.95 

10.91 

0.045 

6 

0.000374 

0.000350 

0.000391 

0.000360 

0.000011 

P(D 

0.8511 

0.8251 

0.8543 

0.8403 

0.0074 

P(2) 

0.6973 

0.6559 

0.7067 

0.6752 

0.0141 

ai 

0.9343 

0.890 

0.961 

0.929 

0.0187 

a2 

-0.0975 

-.136 

-0.078 

-0.106 

0.0191 

Ax  (lyr) 

1.275 

1.217 

1.334 

1.278 

0.036 

A2  (2yr) 

0.402 

0.360 

0.429 

0.401 

0.023 

Since  the  estimated  s.d.    of  the  slope  of  the   (linear)   long  term  trend 
m(t)   =  a+6t     is  0.000011,   it   is  clear  that  the  tilt   in  the  data  is  statistically 
significant.     Note  that  gross  approximations  give  the  s.d. 's  of  the  estimates 
of     p(l)   and  p(2)   in  a  sample  of  size  4380  as     1/(4380)1/2  =  0.015  which, 
given  statistical   fluctuations,   are   in   reasonable  agreement  with  the  estimated 
s.d. 's   of  the  serial   correlations  of  0.0074  and  0.0141.      Note,    however,    that  the 
estimates   of  p(l)   and  p(2)   have  averages   below  the  nominal    values.     This   reflects 
the  known  fact  that  estimates  of  serial   correlations  with  high  values   are  biased; 
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there   is    a  technique  for   removing  this   kind  of   bias   called  jacknifing    (Miller, 
1974)   and  this   could  have  been  used.     Again,   the  estimated  precision  of  the 
estimate  of     A1    is     0.036/1.275     x   100%  =  2.8%,   which   is   suitably   low. 

To  investigate  the  problem  of  sampling  rate,    the  simulations   in  Table   1 
were   run  to  produce  estimates   of     a,    6,   A     and  A     based  on  every  2nd  point   in 
the  generated  series    (Table  8)   and  every  third  point    (Table   9). 

Table  8   (every  2nd  point) 
Parameter         Nominal    Value  Min  Max  Average  s. d. 


a 

10.9 

10.79 

10.95 

10.91 

0.042 

6 

0.000374 

0.000352 

0.000387 

0.000360 

0.000010 

Al 

1.275 

1.215 

1.330 

1.277 

0.035 

A2 

0.402 

0.362 

0.429 

0.400 

0.023 

Table  9  {eyery  3rd  point) 
Parameter         Nominal    Value  Min  Max  Average  s. d. 


a 

10.9 

10.78 

10.95 

10.91 

0.049 

6 

0.000374 

0.000348 

0.000393 

0.000362 

0.000012 

Al 

1.275 

1.213 

1.342 

1.279 

0.037 

A2 

0.402 

0.369 

0.429 

0.403 

0.022 

First   note  that   the  estimated  standard  deviations   in  the  last  columns   in 
Tables   7  and  8  are   virtually     identical.      For   independent  data  the  s.d. 's   in 
Table  8  would  be  approximately     /2     times   those  in  Table   7.     The  results   there- 
fore  indicate  that   day-to-day   correlation   in  the  data   is   so  high    (local   scales 
are  so  long)  that   sampling  every   day  does   not   give  any   improvement   to  the 
sampling  of   long  scale    (one  year  and  up)   phenomena.     Such  phenomena  are 
effectively   summarized  by   the  model   through  s(t)   and  m(t),   exclusive  of  e(t). 
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In  Table  9,  results  are  given  for  estimation  with  sampling  every  third  day. 
The  standard  deviations  are  beginning  to  increase,  but  not  sufficiently  to 
justify  any  finer  sampling. 

10.  Oceanographic  aspects  of  the  data 

10.1  Spring  Transitions 

Abrupt  decreases  in  SST  have  been  observed  in  the  annual  temperature  cycle 
at  Granite  Canyon  during  the  spring  and  appear  to  signal  the  onset  of  upwelling 
along  the  central  California  coast.  These  dramatic  changes  in  temperature  occur 
in  at  least  6  out  of  the  12  years  considered  and  have  been  individually 
identified  and  tabulated  previously  in  Table  1.  From  the  table  it  is  seen  that 
these  transitions  typically  occur  in  March  and  last  for  about  a  week.  Also  the 
associated  temperatures  drop  from  about  13  to  about  IOC.  We  have  examined 
one  spring  transition  in  particular,  the  spring  transition  in  1980  starting  on 
March  11th,  which  is  shown  in  Figure  22.  During  this  transition,  temperatures 
dropped  about  4C  over  a  period  of  7  days.  As  indicated  previously  in  Figure  3, 
the  entire  water  column  experienced  this  transition  down  to  at  least  500  m. 
Alongshore,  currents  measured  at  the  same  location  reversed  from  poleward  to 
equatorward  in  direction  at  each  of  four  depths  between  113  and  510  m  at  the 
same  time  as  the  decreases  in  temperature  occurred.  The  synoptic-scale  winds  at 
this  location  however  did  not  become  upwelling  favorable  until  just  after  this 
event  occurred.  The  latter  observation  suggests  that  this  event  may  have  been 
caused  by  wind  forcing  at  another  location.  During  a  similar  spring  transition 
during  1981,  Brink,  Stuart,  and  Van  Leer  (1983)  also  concluded  that  local  wind 
forcing  was  not  responsible  for  the  transition  they  observed  off  Pt.  Conception 
during  that  year.  Finally  IR  satellite  imagery  acquired  a  few  days  before  and 
immediately  following  the  event  (not  shown)  indicated  a  change  from  isothermal 
surface  conditions  to  conditions  suggestive  of  intense  upwelling  in  the  vicinity 
of  Pt.  Sur. 
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The  relative  times  at  which  the  1980  spring  temperature  transition  arrived 
at  the  four  locations  along  the  California  coast  for  which  we  have  data  are  seen 
in  Figures  4,  5,  and  6.  This  event  is  seen  to  arrive  first  at  Diablo  Canyon, 
second  at  Granite  Canyon,  third  at  the  Farallons,  and  finally  at  Pacific  Grove. 
Phase  speeds  for  poleward  propagating  disturbances  along  the  coasts  of  North  and 
South  America  generally  range  between  about  35  and  280  Km/day.   In  particular, 
Allen  and  Romea  (1980)  indicate  that  phase  speeds  of  the  order  of  175  Km/day 
should  be  expected  for  barotropic  continental  shelf  waves  along  the  South 
American  coast  south  of  about  18°S.  We  take  this  value  as  representative  along 
the  North  American  coast  as  well,  since  the  continental  shelf-slope  regions 
along  both  South  and  North  America  are  generally  similar  (Shepard,  1973).  The 
difference  in  arrival  time  between  Diablo  Canyon  and  Granite  Canyon,  which  are 
about  175Km  apart,  is  approximately  a  day,  a  time  separation  consistent  with  the 
propagation  speed  of  175  Km/day  for  barotropic  continental  shelf  waves  in  this 
region.  However,  the  separation  in  arrival  times  between  Granite  Canyon  and  the 
Farallons  is  considerably  greater  (~  6  days)  even  though  the  distance  between 
these  points  is  again  about  175  Km.  The  time  separation  in  this  case  is  not 
consistent  with  propagation  of  continental  shelf  waves  since  phase  speed  should 
in  fact  increase  due  to  the  slight  increase  in  the  Coriolis  parameter  and  an 
increase  in  width  of  the  continental  shelf  in  this  region.  Consequently  we 
conclude  that  these  results  are  not  necessarily  consistent  with  shelf  wave 
dynamics.  We  also  note  that  the  sea  surface  temperatures  are  only  sampled  once 
a  day.  Thus  our  ability  to  resolve  the  phase  speeds  of  these  waves  is  limited. 

10.2  Annual  cycles 

A  mean  annual  temperature  cycle  has  been  calculated  for  the  Granite  Canyon 
data,  realizing  from  our  previous  statistical  calculations  that  a  mean  value  may 
not  be  a  truly  representative  measure  of  central  tendency  in  our  data.  Points 
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corresponding  to  the  average  temperatures   and  a  smoothed  version  of  that   average 
(with   a  window  of  25  points)   are   both   presented   in  Figure  23.    According  to  List 
and  Koh    (1976),    and  Reid,   Roden  and  Wyllie    (1958),    a  well-defined  annual    cycle 
in  SST  might   not   be  expected  off  Granite  Canyon,    due  to  warming  during  the 
winter  from  the  Davidson  Current   and  cooling  during  the  spring  and  summer  from 
coastal   upwelling.     To  the  contrary,   Figure  23  shows   a  well-defined  annual   cycle 
at  Granite  Canyon,   with   a  range  of  about   3C.   However  the  annual    range  of 
variability   in  temperature  is   reduced  at  Granite  Canyon  compared  to  the  range  of 
variation  usually  found  at   this   latitude.     This   is   shown  in  the  insert   in  Figure 
23  where  deep  ocean  SSTs    (Robinson,    1975)   are   compared  to  the  Granite  Canyon 
SSTs.   The  effect   of  averaging  the  data  has   been  in  part  to   "smear"  the 
occurrence  of  the  spring  transition  over  a  two-month  period  from  about 
mid-February  to  mid-April.     Maximum  temperatures   are  seen  to  occur  during 
mid-October,    about  a  month   later  in  the  year  than  might   be  expected  at   this 
latitude.     Obviously  the  phase   relationships   between  such  factors   as   the  onset 
and  cessation  of  the  Davidson  Current,   and  the  onset  and  cessation  of  coastal 
upwelling  determine   in  part  the  detailed  structure   of  the  annual   temperature 
cycle  off  Pt.   Sur  in  any  given  year. 

Skogsberg   (1936)   originally   specified  three  oceanic  seasons   for  Monterey 
Bay,    an  upwelling  season  from  about  February  to  July,    an  oceanic  season   from 
about  July  to  November,    and  a  Davidson  Current   season  from  about  November  to 
February.     According  to  the  smoothed  average  annual   cycle  presented  in  Figure 
23,   at   least  two   "seasons"  can  be  identified,    an  upwelling  season  and  a 
non-upwel 1 i ng  season.      Based  on  our  data  the  distinction   between  an  oceanic 
season   and  a  Davidson  Current  season   is   not   clear.      We  note  however  that 
Skogsberg  used  spatially  distributed  temperature  data  as  well    as   chemical   and 
biological   observations   to  construct  his   seasonal    description.      Also  his 
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seasonal   model   was   derived  for  Monterey  Bay   and  as   a   result   it  may   not   be   valid 
to  extrapolate  his   model   to  regions   outside  the  Bay. 

10.3  El   Nino  events 

In  addition  to  the  seasonal    variation   in  temperature  at  Granite  Canyon, 
inter-annual    variability   is   also  evident,    as   mentioned  previously   in  Section   3. 
The   large  increases   in  temperature  seen  during  the  winters   of   1972,    1976,   and 
1983     correspond  to  major  El   Nino  warming  events   that  took   place  during  those 
years.   These  events   have  also  been  observed  off  central   California  through 
biological   measurements   and  measurements   of  sea   level    by  Chelton,   Bernal   and 
McGowan    (1982).     The  increase  in  temperature  seen  during  the   latter  part  of   1979 
is   not   associated  with   an  El    Nino     warming  event    (Enfield,    1983). 

Of  particular  interest   is   the  present   major  El   Nino     warming  event.     This 
warming  was   shown  to  extend  at   least   500  Km  off  the  coast  of  central   California 
during  February   1983   (Under  and  Breaker,    1983).     The  effects   of  the  present 
warm  anomaly  are  clearly  evident   in  Figure  24  which   compares   the  3/1/82  to 
3/1/83  time  period    (i.e.,   the  solid  curve)   to  the  12  year  average  data   (i.e., 
the  dashed  curve).    Early   in   1982  warming  off  Granite  Canyon  appeared  to  be 
intermittent.      However,    by   late  1982  warming  at  this    location  had  become 
well-established.     As  of  February   1983,   the  warm  anomaly   at  Granite  Canyon  had 
reached  almost   3C.   A  smoothed  version   of  Figure  24   is   shown  in  Figure  2b  and 
helps   clarify  the  present   departure   in  SST   from  the   12-year  mean  at  this 
location. 

10.4     Interpretation  of  higher  frequency   variations 

From  an  oceanographic  viewpoint,   the  periodogram  shown  previously   in 
Figure   9  is   somewhat   surprising  for   its   lack   of  any  major  spectral   components 
beyond   100  cycles.     For  example,    spectra  of  winds   along  the  California  coast 
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often  show  higher  amplitudes   at  periods   between  3  and   10  days   (Van  Patten, 
1981).     There  are  three  possible  explanations   for  this   apparent   lack   of  con- 
sistency between  wind  and  temperature  spectra.     First,   the  periodic  effects   in 
the  winds  may   be   real,   but  continual   stirring  of  the  near-surface  environment 
where  observations  of  ocean  temperature  are  made,   may  inhibit  these  periodic  in- 
fluences on  SSTs.     Alternatively,   higher  amplitudes   in  the  coastal  wind  spectra 
may  be  caused  by   correlation   in  the  wind  data   rather  than  by  underlying  period- 
icities.    As  such,   a  periodic-like  response  in  SST  to  the  local  winds  should  not 
necessarily  be  expected.     A  third  possibility  of  course  is  that  the  spectrum  of 
the  wind  field  at  this  particular  location  has   no  discrete  components  to  start 
with.     Data  to  check  this  third  possibility  unfortunately  do  not  exist  over  the 
full    12  year  period  but  there  is  evidence  from  other  locations    (Hugus,   1982) 
that  3  to  10  day  periodicities  do  not  show  up  consistently  in  wind  observations. 

Originally  the  possibility  of  aliasing  by  the  semiduirnal   and/or  diurnal 
tides  was   a  concern.     Analysis  of  similar  observations   at  other  locations  along 
the  California  coast  by  List  and  Koh    (1976)   indicated  that  tidal   aliasing   (or  a 
natural   periodicity  corresponding  to  the  spring/neap  tides)  by  the  semidiurnal 
tide  at   14.78  days  could  be  important.     To  the  contrary  we  do  not  find  evidence 
that  the  alias  of  either  tidal   component    (148  cycles  for  the  diurnal  tide  and 
296  cycles  for  the  semidiurnal   tide)  appears   in  the  calculated  periodogram. 

In  examining  the  periodogram  of  second  residuals,   particularly  at  high 
frequencies,    it  is  of   interest  to  briefly  consider  the  likely  influence  of 
turbulence  on  these  results.        Turbulence  production  undoubtedly  occurs   from 
the  breaking  of  incoming  surface  waves  and  perhaps  from  local   drift  currents  as 
well.      In  this   regard  a  portion  of  the  periodogram  of  second  residuals  has  been 
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smoothed  and  replotted  on  log  coordinates  in  Figure  26.  The  uniform  decrease 
in  spectral  amplitude  with  increasing  frequency  is  reminiscent  of  the  spectral 
decay  often  associated  with  the  inertial  subrange  from  homogeneous  turbulence 
theory  (Batchelor,  1953).  Between  200  (~  2  days)  and  2000  cycles  (-22  days), 
the  spectral  slope  is  essentially  constant  at  -y/b,  close  to  the  value  of  -b/3 
predicted  by  classical  turbulence  theory  (Phillips,  1977 ) .  According  to  this 
theory  temperature  is  assumed  to  be  a  conservative  quantity,  an  assumption  which 
could  be  quite  unrealistic  at  the  sea-surface.   In  our  case  time  replaces 
distance  which  is  normally  taken  as  the  independent  variable,  and  if  it  is 
assumed  that  these  coordinate  frames  are  linearly  related,  then  the  spectral 
slopes  observed  here  are  directly  comparable  with  those  usually  derived  in 
oceanic  turbulence  calculations  as  a  function  of  distance.   In  turbulence  models 
where  temperature  is  not  assumed  to  be  conservative,  spectral  slopes 
considerably  greater  than  -b/3  are  predicted  (Panchev,  1980).  Measurements  of 
the  spatial  surface  temperature  spectrum  by  Saunders  (1972)  and  Holladay  and 
O'Brien  (197b)  for  example,  yielded  spectral  slopes  of  -2.2  and  -3  respectively, 
suggesting  that  temperature  should  not  be  treated  as  a  conservative  quantity  at 
least  under  certain  conditions. 

In  conclusion,  our  earlier  analyses  suggest  that  the  (likely)  influence 
of  turbulence  in  our  data  corresponds  to  the  autoregressi ve  component  in  our 
statistical  decomposition  since  the  spectrum  of  the  fourth  residuals  is 
essentially  flat.  A  more  detailed  analysis  of  our  data  would  include  a 
comparison  of  the  spectra  for  various  orders  of  autoregression  with  spectra 
predicted  from  various  theoretical  models  of  ocean  turbulence. 
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